void DFT(double** input, int size)

void FFT(complex* input, int size);
void sFFT(double** input, int size, double epsilon, double delta);

struct complex {
    double v[2];
    complex(){
        v[0] = v[1] = 0.0;
    }
    double mgntd(){
        return v[0][0] * v[0][0] + v[1][1] * v[1][1];
    }
    complex operator*(const complex& b)const{
        complex res;
        res.v[0] = v[0] * b.v[0] - v[1] * b.v[1];
        res.v[1] = v[0] * b.v[1] + v[1] * b.v[0];
        return res;
    }
    complex operator/(const complex& b)const{
        complex res;
        double temp = b.v[0] * b.v[0] + b.v[1] * b.v[1];
        res.v[0] = (v[0] * b.v[0] + v[1] * b.v[1]) / temp;
        res.v[1] = (b.v[0] * v[1] - b.v[1] * v[0]) / temp;
        return res;
    }

}


int divideArray(complex* data, int* indices, int size){
    int iter = -1;
    double pivot = data[size - 1].mgntd();
    for (int i = 0; i < size - 1; ++i){
        if (data[i].mgntd() > pivot){
            ++iter;
            swap(data[i], data[iter]);
            swap(indices[i], indices[iter]);
        }
    }
    swap(data[size - 1], data[iter]);
    swap(indices[size - 1], indices[iter]);
    return iter;
}

void findLargest(complex* data, int* indices, int num, int size){
    
    int p = divideArray(data, indices, size) + 1;
    if (p == num)
        return;
    if (p < num)
        findLargest(data + p, indices + p, num - p, size - p);
    if (p > num)
        findLargest(data, indices, num, p);
            
}

//size = 2^p

int hash(int sigma, int i, int B, int size){
    return sigma * i * B / size;
}

int offset_hash(int sigma, int i, int B, int size){
    int temp =  sigma * i - hash(sigma, i, B, size) * (size / B);

}
int* locationLoop(double** input, int size, complex* z, double* window,
                 int support, int B, int dk,  unsigned sigma, unsigned tau){
    for (int i = 0; i < B; ++i){
        int limit = (support / B) + (support % b != 0) - 1;
        for (int j = 0; j < limit; ++j){
            int index  = (i + j * B) % size;
            z[i].v[0] += window[index] * input[(sigma * index + tau) % size][0];
            z[i].v[1] += window[index] * input[(sigma * index + tau) % size][1];
        }
    }
    FFT(z, B);
    int* indices = new int[B];
    for (int i = 0; i < B; ++i)
        indices[i] = i;
    findLargest(z, indices, dk, B);  

    
    int* result = new int [dk * size / B + 5];
    int iter = 1;
    result[0] = 0; // size    
    for (int i = 0; i < dk; ++i){
        int temp = indices[i] * (size / B) / sigma;
        while (temp > 0 && hash(sigma, temp, B, size) >= indices[i])
            --temp;
        if (hash(sigma, temp, B, size) < indices[i])
            ++temp;
        while (temp < size && hash(sigma, temp, B, size) == indices[i]){
            result[iter++] = temp;
            result[0]++;
        }
    }
    //restore z^
    for (int i = 0; i < B; ++i){
        swap(z[i], z[indices[i]]);
        swap(indices[i], indices[indices[i]]);
    }
    delete[] indices;     
    return result;
}


double* window_function;
int window_support;

complex unity(int n, int sz){
    complex res;
    res.v[0] = cos((2.0 * M_PI / sz) * n);
    res.v[0] = sin((2.0 * M_PI / sz) * n);
    return res;
}

complex* estimationLoop(int* input, int size, fullsz, complex* z, complex* DFTwindow,
                        int B, unsigned sigma, unsigned tau){
    complex* out = new complex[size];
    for (int i = 1; i <= size; ++i){
        int j = i - 1;
    
        out[j] = z[hash(sigma, input[i], B, fullsz)] * unity(tau*input[i], fullsz) / DFTWindow[offset_hash(sigma, input[i], B, fullsz)];
    }    
    return out;
}
set<int,  sFFT(double** input, int size, int sparsity, double epsilon, double delta){
    int B = 10;
    double d = 1.0 / epsilon * 1.0;
    int L = 1; // O(log N);
    unsigned* sigma = new unsigned[L];
    unsigned* tau = new unsigned[L];
    complex z** = new complex*[L];    
    for (int i = 0; i < L; ++i){
        z[i] = new complex[B];
        sigma[i] = (rand() % (size / 2)) * 2 + 1;
        tau[i] = rand() % size;  
    }
    int** indices_set = new int*[L];
    for (int i = 0; i < L; ++i)
        indices_set[i] = locationLoop(input, size, z[i], window_function, window_support, B, d * sparsity, sigma[i], tau[i]);
    map<int, int> indices_frequency;
    for (int i = 0; i < L; ++i)
        for (int j = 1; j <= indices_set[i][0];  ++j){
            if (indices_frequency.find(indices_set[i][j]) == indices_frequency.end())
                indices_frequency[indices_set[i][j]] = 1;
            else
                indices_frequency[indices_set[i][j]]++;
        }
    int cnum = indices_frequency.size();
    int* big_coeffs = new int[cnum];
    int iter = 0;
    for (auto i = indices_frequency.begin(); i != indices_frequency.end(); ++i)
        if (i->second >= L / 2)
            big_coeffs[iter++] = i->first;

    complex** estimates = new complex*[cnum];
    for (int i = 0; i < cnum; ++i)
        estimates[i] = estimationLoop(big_coeffs, cnum, size, z[i], DFTwindow, B, sigma[i], tau[i]);


    double* candidatesRe = new double[L]
    double* candidatesIm = new double[L]
    unordered_set<int, complex> res;
    for (int i = 0; i < cnum; ++i){
        for (int j = 0; j < L; ++j){
            candidatesRe[j] = estimates[j][i].v[0];
            candidatesIm[j] = estimates[j][i].v[1];
        }
        sort(candidatesRe, candidatesRe + L);//Can use linear algorithm instead, but that would be log log n speed up, so who cares?
        sort(candidatesIm, candidatesIm + L);
        complex t;
        t.v[0] = candidatesRe[L / 2];
        t.v[1] = candidatesIm[L / 2];
        res[big_coeffs[i]] = t;
    }

    delete[] candidatesRe;
    delete[] candidatesIm;
    delete[] tau;
    delete[] sigma;
    for (int i = 0; i < L; ++i){
        delete[] z[i];
        delete[] indices_set[i];
    }
    delete[] z;
    delete[] indices_set;
    delete[] big_coeffs;
    for (int i = 0; i < cnum; ++i)
        delete[] estimates[i];
    delete[] estimates;
    return res;
}